{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 古典线性回归"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 导入样本数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>日期</th>\n",
       "      <th>hs300</th>\n",
       "      <th>sz</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2019-11-22</td>\n",
       "      <td>3849.9948</td>\n",
       "      <td>2885.2884</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2019-11-21</td>\n",
       "      <td>3889.5980</td>\n",
       "      <td>2903.6379</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2019-11-20</td>\n",
       "      <td>3907.8641</td>\n",
       "      <td>2911.0534</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2019-11-19</td>\n",
       "      <td>3947.0392</td>\n",
       "      <td>2933.9908</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2019-11-18</td>\n",
       "      <td>3907.9291</td>\n",
       "      <td>2909.2002</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "          日期      hs300         sz\n",
       "0 2019-11-22  3849.9948  2885.2884\n",
       "1 2019-11-21  3889.5980  2903.6379\n",
       "2 2019-11-20  3907.8641  2911.0534\n",
       "3 2019-11-19  3947.0392  2933.9908\n",
       "4 2019-11-18  3907.9291  2909.2002"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "\n",
    "df = pd.read_excel('./数据/上证指数与沪深300.xlsx')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 原理讲解"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 古典线性回归模型的假定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 假定1：线性假定"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "总体模型为:$y_i=\\beta_1x_{i1}+\\beta_2x_{i2}+\\cdots+\\beta_Kx_{iK}+\\varepsilon_i (i=1,\\cdots,n)\\tag{1}$            \n",
    "\n",
    "其中 $n$ 为样本容量，解释变量 $x_{ik}$ 的第一个下标表示第$i$个“观测值”，而第二个下标则表示第 $k$ 个解释变量$(k=1,\\cdots,K)$，共有 $K$个 解释变量。如果有常数项，则通常令第一个解释变量为单位向量，即 $x_{i1}=1$\n",
    "\n",
    "为了更简洁地表达，下面引入矩阵符号。把方程（1）的所有解释变量和参数都写成向量，记第$i$个观测数据为$x_i\\equiv\\left(x_{i 1} ,x_{i_{2}} \\cdots x_{i K}\\right)^{\\prime}$，$\\beta \\equiv\\left(\\beta_{1}, \\beta_{2} \\cdots \\beta_{K}\\right)^{\\prime}$，则方程（1）为：\n",
    "\n",
    "$y_{i}=x_{i}^{\\prime} \\boldsymbol{\\beta}+\\varepsilon_{i} \\quad(i=1, \\cdots, n)\\tag{2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 假定2：严格外生性"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\mathrm{E}\\left(\\varepsilon_{i} | X\\right)=\\mathrm{E}\\left(\\boldsymbol{\\varepsilon}_{i} | x_{1}, \\cdots, x_{n}\\right)=0 \\quad(i=1, \\cdots, n)\\tag{3}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 假定3：不存在“严格多重共线性”（strict multicolinearity），即数据矩阵 $X$ 满列秩"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 假定4：球型扰动项（spherical disturbance），即扰动项满足“同方差”、“无自相关”的性质"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\operatorname{Var}(\\boldsymbol{\\varepsilon} | X)=\\mathrm{E}\\left(\\boldsymbol{\\varepsilon} \\boldsymbol{\\varepsilon}^{\\prime} | X\\right)=\\sigma^{2} \\boldsymbol{I}_{n}=\\left(\\begin{array}{ccc}\\sigma^{2} & & 0 \\\\ & \\ddots & \\\\ 0 & & \\sigma^{2}\\end{array}\\right)\\tag{4}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 最小二乘法\n",
    "假定待估计方程为：$hs300 = c+sz$。其中 c 为常数项"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### OLS估计量b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$b \\equiv\\left(X^{\\prime} X\\right)^{-1} X^{\\prime} y\\tag{5}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OLS估计值为：\n",
      " [[-124.69031687]\n",
      " [   1.29305435]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "n = df.shape[0]   # 样本容量\n",
    "beta = np.array(df['sz']).reshape(n,1)\n",
    "c = np.ones((n,1))   # 常数项\n",
    "X = np.hstack((c,beta))   # hstack()在行上合并，vstack()在列上合并\n",
    "y = np.array(df['hs300']).reshape(n,1)\n",
    "\n",
    "b = np.linalg.inv(X.T @ X) @ X.T @ y\n",
    "print('OLS估计值为：\\n',b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 残差 $e$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$e \\equiv\\left(\\begin{array}{llll}\n",
    "e_{1},e_{2},\\cdots,e_{n}\n",
    "\\end{array}\\right)=y-X \\widetilde{\\beta}\\tag{6}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "e = y - X @ b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 扰动项方差 $s^{2}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "对于扰动项方差$\\sigma^{2}=\\operatorname{Var}\\left(\\varepsilon_{i}\\right)$，由于总体扰动项 $\\varepsilon$ 不可观测，而样本残差 $e$ 可以近似地看成是 $\\varepsilon$ 的实现值，故使用以下统计量作为对方差 $\\sigma^{2}$ 的估计：\n",
    "$$s^{2} \\equiv \\frac{1}{n-K} \\sum_{i=1}^{n} e_{i}^{2}\\tag{7}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "平方和： [[5453855.41555515]]\n",
      "扰动项方差 [[11907.98125667]]\n",
      "扰动项标准差 109.12369704454953\n"
     ]
    }
   ],
   "source": [
    "K = X.ndim\n",
    "SSE = e.T @ e\n",
    "s2 = SSE/(n-K) \n",
    "\n",
    "import math \n",
    "s = math.sqrt(s2)\n",
    "\n",
    "print('平方和：', SSE)\n",
    "print('扰动项方差', s2)\n",
    "print('扰动项标准差', s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 估计量 b 的方差-协方差矩阵"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\\operatorname{Var}(b | X)=\\boldsymbol{\\sigma}^{2}\\left(X^{\\prime} X\\right)^{-1}\\tag{8}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "协方差矩阵：\n",
      " [[ 3.77128774e+03 -1.27819004e+00]\n",
      " [-1.27819004e+00  4.36206925e-04]]\n"
     ]
    }
   ],
   "source": [
    "Varb = s2 * np.linalg.inv(X.T @ X)\n",
    "\n",
    "print('协方差矩阵：\\n', Varb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 置信区间"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "由于 $\\frac{b_{k}-\\beta_{k}}{\\mathrm{SE}\\left(b_{k}\\right)} \\sim t(n-K)$，根据 $t_{\\alpha/2}$ 得：\n",
    "$$\\mathrm{P}\\left\\{-t_{\\alpha / 2}<\\frac{b_{k}-\\beta_{k}}{\\mathrm{SE}\\left(b_{k}\\right)}<t_{\\alpha / 2}\\right\\}=1-\\alpha\\tag{9}$$\n",
    "\n",
    "$$P\\left\\{b_{k}-t_{\\alpha / 2} \\operatorname{SE}\\left(b_{k}\\right)<\\boldsymbol{\\beta}_{k}<b_{k}+t_{\\alpha / 2} \\operatorname{SE}\\left(b_{k}\\right)\\right\\}=1-\\alpha\\tag{10}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "95% 置信区间：\n",
      " [[-245.37220852   -4.00842522]\n",
      " [   1.25201092    1.33409777]]\n"
     ]
    }
   ],
   "source": [
    "from scipy. stats import t\n",
    "\n",
    "alpha = 0.05   # 置信度\n",
    "nu = max(0,n-K)   # 自由度\n",
    "tval = t.ppf(1-alpha/2,nu)   # 逆函数值\n",
    "SE_b = np.sqrt(np.diag(Varb)).reshape(K,1)\n",
    "bint = np.hstack((b-tval*SE_b,b+tval*SE_b))\n",
    "print('95% 置信区间：\\n', bint)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### t检验"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$t_{k} \\equiv \\frac{b_{k}-\\bar{\\beta}_{k}}{\\mathrm{SE}\\left(b_{k}\\right)} \\equiv \\frac{b_{k}-\\bar{\\beta}_{k}}{\\sqrt{s^{2}\\left(X^{\\prime} X\\right)_{k k}^{-1}}} \\sim t(n-K)\\tag{11}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div align=center><img src=\"https://lei-picture.oss-cn-beijing.aliyuncs.com/img/20200423133127.png\" width=\"450\"></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "t检验为：\n",
      " [[-2.0304294 ]\n",
      " [61.91138223]]\n",
      "\n",
      "\n",
      "p值为：\n",
      " [[0.0428902]\n",
      " [0.       ]]\n"
     ]
    }
   ],
   "source": [
    "t_stat = b/SE_b\n",
    "t_p = 2*(1-t.cdf(abs(t_stat),n-K))\n",
    "\n",
    "print('t检验为：\\n', t_stat)\n",
    "print('\\n')\n",
    "print('p值为：\\n', t_p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 两类错误"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "根据样本信息对总体进行推断，有可能犯错误。特别地，在进行假设检验时，可能犯两类性质不同的错误。\n",
    "\n",
    "**第Ⅰ类错误**：虽然原假设为真，但却根据观测数据做出了拒绝原假设的错误判断，即“弃真”。\n",
    "\n",
    "**第Ⅱ类错误**：虽然原假设为假（替代假设为真），但却根据观测数据做出了接受原假设的错误判断，即“存伪”。\n",
    "\n",
    "由于在进行假设检验时，通常知道第Ⅰ类错误的发生概率，而不知道第Ⅱ类错误的发生概率。因此，如果拒绝原假设，可以比较理直气壮，因为知道犯错误的概率就是显著性水平（比如5%）；另一方面，如果接受原假设，则比较没有把握，因为我们通常并不知道第Ⅱ类错误的发生概率（可能很高）。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 使用 statsmodels 库实现"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"simpletable\">\n",
       "<caption>OLS Regression Results</caption>\n",
       "<tr>\n",
       "  <th>Dep. Variable:</th>            <td>y</td>        <th>  R-squared:         </th> <td>   0.893</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Model:</th>                   <td>OLS</td>       <th>  Adj. R-squared:    </th> <td>   0.893</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Method:</th>             <td>Least Squares</td>  <th>  F-statistic:       </th> <td>   3833.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Date:</th>             <td>Thu, 23 Apr 2020</td> <th>  Prob (F-statistic):</th> <td>1.20e-224</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Time:</th>                 <td>14:24:33</td>     <th>  Log-Likelihood:    </th> <td> -2810.3</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>No. Observations:</th>      <td>   460</td>      <th>  AIC:               </th> <td>   5625.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Residuals:</th>          <td>   458</td>      <th>  BIC:               </th> <td>   5633.</td> \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Df Model:</th>              <td>     1</td>      <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Covariance Type:</th>      <td>nonrobust</td>    <th>                     </th>     <td> </td>    \n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "    <td></td>       <th>coef</th>     <th>std err</th>      <th>t</th>      <th>P>|t|</th>  <th>[0.025</th>    <th>0.975]</th>  \n",
       "</tr>\n",
       "<tr>\n",
       "  <th>const</th> <td> -124.6903</td> <td>   61.411</td> <td>   -2.030</td> <td> 0.043</td> <td> -245.372</td> <td>   -4.008</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>x1</th>    <td>    1.2931</td> <td>    0.021</td> <td>   61.911</td> <td> 0.000</td> <td>    1.252</td> <td>    1.334</td>\n",
       "</tr>\n",
       "</table>\n",
       "<table class=\"simpletable\">\n",
       "<tr>\n",
       "  <th>Omnibus:</th>       <td>61.627</td> <th>  Durbin-Watson:     </th> <td>   0.010</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Prob(Omnibus):</th> <td> 0.000</td> <th>  Jarque-Bera (JB):  </th> <td>  83.387</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Skew:</th>          <td> 1.031</td> <th>  Prob(JB):          </th> <td>7.81e-19</td>\n",
       "</tr>\n",
       "<tr>\n",
       "  <th>Kurtosis:</th>      <td> 2.692</td> <th>  Cond. No.          </th> <td>3.55e+04</td>\n",
       "</tr>\n",
       "</table><br/><br/>Warnings:<br/>[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.<br/>[2] The condition number is large, 3.55e+04. This might indicate that there are<br/>strong multicollinearity or other numerical problems."
      ],
      "text/plain": [
       "<class 'statsmodels.iolib.summary.Summary'>\n",
       "\"\"\"\n",
       "                            OLS Regression Results                            \n",
       "==============================================================================\n",
       "Dep. Variable:                      y   R-squared:                       0.893\n",
       "Model:                            OLS   Adj. R-squared:                  0.893\n",
       "Method:                 Least Squares   F-statistic:                     3833.\n",
       "Date:                Thu, 23 Apr 2020   Prob (F-statistic):          1.20e-224\n",
       "Time:                        14:24:33   Log-Likelihood:                -2810.3\n",
       "No. Observations:                 460   AIC:                             5625.\n",
       "Df Residuals:                     458   BIC:                             5633.\n",
       "Df Model:                           1                                         \n",
       "Covariance Type:            nonrobust                                         \n",
       "==============================================================================\n",
       "                 coef    std err          t      P>|t|      [0.025      0.975]\n",
       "------------------------------------------------------------------------------\n",
       "const       -124.6903     61.411     -2.030      0.043    -245.372      -4.008\n",
       "x1             1.2931      0.021     61.911      0.000       1.252       1.334\n",
       "==============================================================================\n",
       "Omnibus:                       61.627   Durbin-Watson:                   0.010\n",
       "Prob(Omnibus):                  0.000   Jarque-Bera (JB):               83.387\n",
       "Skew:                           1.031   Prob(JB):                     7.81e-19\n",
       "Kurtosis:                       2.692   Cond. No.                     3.55e+04\n",
       "==============================================================================\n",
       "\n",
       "Warnings:\n",
       "[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.\n",
       "[2] The condition number is large, 3.55e+04. This might indicate that there are\n",
       "strong multicollinearity or other numerical problems.\n",
       "\"\"\""
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import statsmodels.api as sm\n",
    "\n",
    "# sm.add_constant(data, prepend=False)\n",
    "mod = sm.OLS(y, X)\n",
    "res = mod.fit()\n",
    "res.summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## matlab实现"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```matlab\n",
    "function [B,resid,siga2,bint,cov_matrix,t,t_p] = OLS_regress(Y,X)\n",
    "\n",
    "%输入变量:\n",
    "%Y - 被解释变量\n",
    "%X - 解释变量\n",
    "\n",
    "%输出变量:\n",
    "% B - 待估计参数beta\n",
    "% resid - 残差\n",
    "% siga2 - 残差方差\n",
    "% bint - 95%置信区间序列\n",
    "% cov_matrix - 协方差矩阵\n",
    "\n",
    "% 1.求OLS估计量B\n",
    "B = inv(X'*X)*X'*Y;\n",
    "\n",
    "% 2.计算协方差矩阵\n",
    "resid = Y - X*B;  %残差\n",
    "[n,K] = size(X);  \n",
    "siga2 = sum(resid.^2)/(n-K);\n",
    "cov_matrix = siga2*inv(X'*X);\n",
    "\n",
    "% 3.t检验\n",
    "t = B./sqrt(diag(cov_matrix));\n",
    "t_p = 2*(1-tcdf(abs(t),n-K));\n",
    "\n",
    "% 4.计算95%置信区间\n",
    "alpha = 0.05;   %置信度\n",
    "nu = max(0,n-K);  %自由度\n",
    "tval = tinv(1-alpha/2,nu);\n",
    "se = sqrt(diag(cov_matrix));\n",
    "bint = [B-tval*se, B+tval*se];\n",
    "```"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
